The thermal Casimir effect in lipid bilayer tubules 
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We calculate the thermal Casimir effect for a dielectric tube of radius R and thickness S formed 
from a membrane in water. The method uses a field-theoretic approach in the grand canonical 
ensemble. The leading contribution to the Casimir free energy behaves as —kBTLnc/R giving rise 
to an attractive force which tends to contract the tube. We find that kc ~ 0.3 for the case of typical 
lipid membrane t-tubules. We conclude that except in the case of a very soft membrane this force 
is insufficient to stabilize such tubes against the bending stress which tends to increase the radius. 



Lipid bilayers in water exhibit a huge variety of ge- 
ometries and in cell biology even more varied structures 
arise. In order to understand where biological mecha- 
nisms, such as molecular motors and cytoskeletal struc- 
tures, are determinant in the stability of biological struc- 
tures, one must first understand the role of the basic 
physical interactions in pure model membrane systems. 
There has been much study of lipid bilayer shape us- 
ing standard continuum mechanics Q, @ . This basic ap- 
proach is also complemented by more microscopic studies 
based on lipid-lipid interaction models 0, U, this ap- 
proach is of course ultimately necessary to fully under- 
stand the physics of bilayers. 

In certain muscle cells, structures known as t-tubules 
are found. These are basically cylindrical tubes whose 
surface is composed of a lipid bilayer. Similar structures 
may also be mechanically drawn off bilayer vesicles. The 
stability of these tubular structures requires an expla- 
nation. The basic continuum theory (where the native 
curvature is zero) Q, predicts that the free energy of a 
tube of length L and radius R is 



F B (L,R) 



R 



(1) 



where the above expression is strictly speaking the excess 
free energy with respect to a flat membrane of the same 
area A — 2ttRL and the subscript B refers to mechanical 
bending. Various experimental and theoretical estimates 
for the bending rigidity Kb can be found in the literature 
and they lie between 1 and 30 (0i0])- The values of k,b 
depend of course on the composition of the bilayer and 
on the experimental protocol used to measure it. One 
crucial clement in both theoretical and experimental de- 
terminations of kb is whether the tube is attached to a 
reservoir of lipid or not, i.e. whether the statistical en- 
semble is grand canonical or grand canonical. Clearly if 
there is no reservoir then any increase in the surface area 
of the tube will lead to a less dense lipid surface concen- 
tration, in this case water may be able to become in con- 
tact with the internal layer composed of the hydrophobic 
heads and so lead to a significant increase in free energy. 
If upon changing the area of the tube lipids can flow into 
the tube wall to maintain the local optimal packing, then 
the free energy cost will be substantially different. This 



bending free energy is positive and hence in the absence 
of external forces or extra constraints it tends to make a 
tube structure expand. We note that when lipid tubes 
are drawn from a vesicle the mechanically applied ten- 
sion can of course overcome this free energy barrier. A 
natural question, motivated by the fact that such struc- 
tures occur in cells, is whether other physical mechanisms 
could lead to their formation and explain their stability. 
A possible explanation put forward is that that electro- 
static effects involving surface charges and ions (salt) in 
the surrounding medium could play a role 0. 0- Cer- 
tain experiments however showed a relative insensitivity 
of the equilibrium tube structure to the concentration of 
salt |5j- There are, however, systems with highly charged 
head groups where the salt concentration does appear 
important in determining the stability of the tubules • 
Another explanation has been put forward in terms of the 
geometry of the lipid, notably the tail having a structure 
such that there is a preferred orientation of the tails next 
to each other, givi ng rise to a chirality which stabilizes 
the tubes |E M E3 ■ 

In this paper we investigate the possibility that a static 
van der Waals or thermal Casimir force could provide an 
attractive force across the tube leading to its stabiliza- 
tion. We adopt a continuum model where the lipid bi- 
layer is modeled as a layer of thickness 6 « 5 — 10(nm) and 
of dielectric constant cm ~ 2eo. The surrounding water 
is also treated as a dielectric continuum of dielectric con- 
stant ew w 80eo. We shall also adopt a model where the 
lipid tube is fixed at each end to a flat lipid reservoir and 
so work in the grand canonical ensemble. The behavior 
of systems composed of layers of varying dielectric con- 
stants was first studied by Lifshitz and coworkers [I3.[T^ . 
The formalism developed is an elegant way of taking into 
account van der Waals forces in a continuum theory. Two 
types of van der Waals forces are accounted for in the- 
ses theories, firstly zero frequency van der Waals forces 
whose nature is classical and secondly the frequency de- 
pendent ones due to temporal dipole fluctuations. In 
terms of thermal field theory the former correspond to 
modes with zero Matsubara frequency and the latter to 
the modes of non-zero frequencies. These latter terms re- 
quire information about the frequency dependence of the 
dielectric constants, where as the former only requires 
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the static dielectric permittivity. In this paper we will 
calculate the contribution of the zero-frequency mode, 
alternatively known as the thermal Casimir effect. The 
quantum Casimir effect corresponding to the modifica- 
tion of the ground state energy of the electromagnetic 
field has been intensively studied in the case of idealized 
boundary conditions in a variety of geometries including 
spheres and cylinders |15| . The thermal Casimir effect 
investigated here has a similar mathematical structure 
though the corresponding effective spatial dimension is 
one less. The temperature dependence of the full Casimir 
effect in a simplified model of a solid dielectric cylinder 
(and sphere) has been recently examined using a heat 
kernel coefficient expansion 16] . In our analysis of the 
diffuse limits we make use of summation theorems for 
Bessel functions which were introduced for the study of 
the Casimir energy for cylinders with light-velocity con- 
serving boundary conditions [TtI ]. 

Wc find that the thermal Casimir effect gives a contri- 
bution to the excess free energy of the tube relative to 
that of the flat plane of 



F C (L,R) 
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[3 log(L<5) + 0.02954] + A i B(A 2 



(2) 



(3) 



with A = (eyy— (-m)/{^w+^m) and where numerically we 
find that B(A 2 ) is only slowly varying with B(0) ~ 1/8. 
Here L « 7r/a is the large wave-number, ultra-violet, 
cut-off which is phenomenologically expressed in terms 
of microscopic cut-off a corresponding to the molecular 
size below which the continuum picture of the dielectric 
medium breaks down. It is not a priori clear whether 
the cut off a should be associated with the water or the 
lipid or indeed both; we discuss this point later. In this 
expression the coefficient of A 2 is dominated by the log- 
arithmic term, the constant contribution being typically 
only of order 1%, and the (estimated) 0(A 4 ) term is of 
comparable size for large enough A. 

We note that the sign of Fc is negative and that it has 
the same functional form as the bending free energy Fb , 
meaning that the thermal Casimir force tends to collapse 
the tube and so helps to stabilize the system against the 
bending forces. We shall show later that with reason- 
able physical parameters kc ~ 0.3 — 0.5. We conclude 
that it is unlikely that the Casimir attraction is able to 
overcome the repulsion due bending that is predicted by 
current theories and measured by experiment. However 
the results of this letter are important for several reasons: 

• We show that the thermal Casimir effect tends to con- 
tract the tube structure. 

• The presence of the microscopic cut off in kc shows that 
the physics is ultimately dominated by the short scale 
or ultra-violet physics. This means that weak electrolyte 



concentrations will have little effect on the system as ver- 
ified by experiments |5j given that there are no strong 
surface charges. 

• Further attractive interactions will be generated by the 
non-zero frequency Matsubara modes. 

• At a technical level use the Pauli van Vleck formula 
to evaluate the arising single body path integrals. This 
technique is direct and well suited to layered systems. 

In what follows we shall sketch the calculation of the 
thermal Casimir free energy in the absence of electrolyte. 
The full calculation is rather lengthy and it will be pre- 
sented in the general case with electrolyte in a more de- 
tailed longer paper 0] . 

The partition function for the zero frequency (static) 
fluctuations of electrostatic field in a medium of varying 
dielectric constant e(x) is given by [l^ ] 



Z = J d[4>] exp ^— ^ J e?xe(x 



)(V0) : 



(4) 



In the model we shall consider the function e(x) = e(r) 
where r is the radial coordinate from the center of the 
tube. The total length of the tube is L and the co- 
ordinate along the length of the tube is denoted by 
z G [-L/2, L/2]. At the extremities of the tube are two 
flat bilayers which act as a reservoir for the material from 
which the tube is made. We consider the limit L > 1 so 
as to be able to neglect edge effects where the tube meets 
the reservoir. The interior of the lipid tube is taken to 
be at Ri = R — 8/2 and the exterior at R2 = R + 6/2 In 
terms of the coordinate r we thus have 



e(r) = e\y r < Ri and r > R2 
e(r) = e M Ri < r < R 2 . 



(5) 



Using the radial symmetry of the problem we express 
the field in cylindrical coordinates as 

00 

(j)(r,z,6)= X rii k(r)exp(ikz)exp(in9) (6) 

n— — 00 k 

The functional integral now factorizes into a product of 
non-interacting path integrals. Taking the limit of L 
large we find 

/elk 
— In (Z nyk (R, e w , cm, $)) 

(7) 



where 



Z nik (R,e w , e M ,S) = J d[X]exp (—^ J 



e(r) 



dX\ 2 I ,r 

-*) + 



k X 



(8) 



As we are only interested in the dependence of the 
free energy as a function of R we may work with the 
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free energy normalized with respect to that of pure water 
with no tube, in this model it is equivalent to subtracting 
off the free energy of a system where ew = £m- This 
regularized free energy is thus 



F(R, e w ,e M ,6) 



k B TL 



dk / Z nt k(R, £w, cm, S) 
2ix \Z rhk (R,e w ,ew,S) 



(9) 



In order to evaluate the individual kernels we make 
use of the Pauli van Vleck formula which is exact for 
quadratic actions. If S[t\, t 2 , X] is an action quadratic in 
X, as is the case here, the Pauli van Vleck formula gives 



K (ti,t 2 ;X,Y) 



X(t 2 )=Y 



X(*i)=X 



d[X] exp(-S[t 1 ,t 2 ,X}) 



1 d 2 S[t 1 ,t 2 ,X cl ] 
' 2ir dXdY 



expi-Sfat^Xd]), (10) 



where X c i is the classical path minimizing the action 
S[t 1 ,t 2 ,X] given by SS[h,t 2 , X]/5X = and with the 
condition that the end points are fixed at X and Y. If 
the action S is quadratic then S[ti, t 2 , X c j\) is a quadratic 
form in X and Y and it only remains to evaluate the re- 
sulting (ordinary) Gaussian integrals after obtainin g th e 
necessary expressions for S[ti, t 2 , X c i\), one obtains [lg 



F(R, ew, £Mj o) 
Lk B T 



^-g(AR u A) + ^-g(AR 2 ,-A) 

ill tt 2 

h(R u R 2 ,A,A) + m(A,A) , (11) 



where A = ir/a is the ultra-violet cut-off corresponding 
to the length scale a. In Eq. Ijlll) we have 



g(x,A) = ^-Y / [ du\n[l + Au(I n (u)K n {u))'} 

h(R 1 ,R 2 ,A,A) = l f — y>[l+ 

{AA 2 k 2 R 1 R 2 F n (kR 1 )I n (kR 1 )K , n (kR 2 )K n (kR 2 ) / 
((1 + Afci? 1 (/„(fci? 1 )A'„(fci? 1 ))') x 
(1 - AkR 2 (I n (kR 2 )K n (kR 2 ))')}] 



m(A 



(12) 



where /„ and A' n are the modified Bessel functions |20j. 
The contribution Lk B Tg(AR, A)/R is the free energy of 
an isolated cylinder of length L, radius R and dielec- 
tric constant cm in a medium of dielectric constant ew- 
Thus the first two terms in Eq. I|ll|l are the respective 
separate contributions of the inner and outer cylindrical 
regions that form the layer of thickness S = R 2 — R\\ the 
term Lk B Th{R\, R 2 , A, A) is the contribution from the 
interaction between the cylinders. The function g(x, A) 



diverges as x — > oo and so this term in the free energy 
must be regulated by taking a non-zero cut-off a. Viewed 
as a Taylor expansion in A we find that the 0(A) term 
of g is independent of r and thus in the free energy the 
contributions proportional to A cancel. This to be ex- 
pected on physical grounds since by examining the limit 
of a diffuse system one can see that any term propor- 
tional to A must be a self energy term [18| . The term of 
order A 2 can be evaluated to be 



9{x,A) 



256 



A 2 [61og(a;) + 30 log 2 + 67 - 11] + 0(A 4 ) 

(13) 



This term is the leading term in the diffuse limit and 
the derivation uses Bessel function summation theorems 
exploited in |17| . The same form for the high temperature 
expansion for a solid dielectric cylinder was obtained in 



expa 

m 



The function h is finite in this same limit and the term 
to, though divergent, is independent of R and cancels 
when the free energy of the bulk in the reservoir is sub- 
tracted. The 0(A 2 ) contribution to h(Ri, R 2 , A, A) is 
given by 



h(R - 5/2, R + 8/2, A, A) = 
3 A 2 1 - y 2 f°° J y 2 z i - 1 



64 R y 2 



dz 



(l + Z 2 2 / 2 ) 1 /2(l + z 2)5/2 



,(14) 



where y = 5/2R. 

Taking account of the bulk reservoir and imposing 
overall conservation of surface area, the relevant free en- 
ergy is that of the tube less the free energy of a flat 
membrane of the same area 

F C {R, ew,£M,S) = F(R, e W) e M , 5) - 2irRLF f , (15) 

where Ff is the free energy per unit area of a flat mem- 
brane, which is be given by 



7T~T „ um ^F(R,e w ,e M ,S), 
k B T ^ A 2m 



167T<5 2 ^— ' TO 3 
m— 1 



(16) 



and can be computed directly in a planar geometry by a 
variety of methods |14| , including the path integral tech- 
nique used here 0] . The contribution to F in this limit 
arises only from the functions h and to in Eq. i|ll|) and 
the latter cancels identically. After subtraction, the dom- 
inant contribution, hc(R, y, A, A), from h is the leading 
1/R term, and we find to 0(A 2 ): 



h c (R,y,A,A) 



j^A 2 
64 ~R 



(log(y) + 21og2-l/2) . (17) 



Putting our results together we find that kq, as defined 
by Eq. J5J, is given by 



KC 



A 2 
64 
A 4 B(A 2 







3 log 


(?) 







61og2 + 37-4 



(18) 



4 



A 

A 


o/a 


0(A 2 ) coeff. of 1/R 
from Eq. ltT£l 


Coeff. of 1/R 
from numerics 


n / A 2 \ 

B(A ) 


78/82 


10 3 


-0.342 


-0.443 


0.123 


78/82 


10 2 


-0.244 


-0.346 


0.123 


0.6 


10 3 


-0.1361 


-0.1520 


0.123 


0.6 


10 2 


-0.0972 


-0.0162 


0.123 


0.2 


10 3 


-0.0151 


-0.0162 




0.6 


10 3 


-0.0038 


-0.0040 





TABLE I: For various values of A and S/a we compare the 
prediction of Eq. I|18[l with numerical integration and deduce 
a numerical value for B(A 2 ). Owing to small systematic er- 
rors in the numerical calculation of the Bessel functions there 
is a negligible discrepancy for very small A but B(A 2 ) is seen 
to be a constant function from evaluations at larger A. We 
see that the result for Fc from Eq. p8p is in very good agree- 
ment with the full calculation. Various values of 8 and a were 
used but typically 5 = 1 — 10(nm) 



where the constant inside the bracket is evaluated to be 
0.02954 . . .. An important point here is the exact cancel- 
lation of the ln(i?) terms coming from g and h, giving a 
leading order behavior of Fc ~ 1/R. The function B (A 2 ) 
receives contribution from both the g and h terms in Eq. 
(|12[) with B(0) ^ 0. Note that there are no odd terms 
in A in the leading 1/R behavior of Fc as, to leading 
order, one may set S/R = in the leading order behavior 
of g and in the denominator of the second term in the 
logarithm of the integral defining h. This is a consistent 
parametrization whilst S S> a. The limit 5 — > must 
be taken carefully and when S < a the separation of F 
in Eq. 112fl into contributions from functions g and h is 
not useful since h develops the compensating ultra-violet 
divergence to that in g and we find lim^o Fc = 0, as 
expected; in essence, the larger of S and a acts as the 
ultra-violet cut-off on the integral defining h. We note 
that the dependence of kc on the membrane thickness 
5 is considerably different to that predicted for the me- 
chanical bending energy, which is 0, El Kb ~ S 2 and is 
verified experimentally [21j if an offset of the measured 
lipid layer thickness 8 is used. Of course Kb as measured 
experimentally should be k^ + kc an d the form predicted 



here for kc is compatible with the qualitative behavior 
seen in measurements plf . 

In Table U for various values of A and S/a, we com- 
pare the prediction of Eq. (|18|) with the result of numer- 
ical integration and deduce a numerical value for B(A 2 ). 
Owing to small systematic errors in the numerical calcu- 
lation of the Bessel functions there is a tiny discrepancy 
for very small A but B(A 2 ) is seen to be a constant func- 
tion from evaluations at larger A and we see that B(0) 
is plausibly 1/8. 

The physical value of the ultra-violet cut-off length can 
only be determined phenomenologically. This is because 
the model is an effective field theory in which the dy- 
namics of the molecular electric dipoles is described by 
the dielectric constant which is a static long-range pa- 
rameter. The field modes with large-A; and n probe the 
static short distance properties of the model and so a 
more refined field theory is needed for these scales. It is 
unclear whether the molecular nature of the lipid has an 
effect on the ultra-violet cutoff but it would seem most 
likely that the effective value of ew at short scales (i.e. 
the microscopic details of water) arc dominant in this 
calculation. 

From Table [I] we see that for a lipid bilayer tube in wa- 
ter with 6 = 10(nm) and a = O.l(nrn) we find kc — 0.346. 
If we were to include the contributions from the modes 
with non-zero Matsubara frequencies, a calculation in 
progress, we can expected at most a factor of two or so 
enhancement based on past experience of similar calcula- 
tions and so kc ~ 1 is a likely largest value. These 
magnitudes are at the lowest end of those for Kb for 
known lipid bilayers in water (0i 0)- However, Wiirger 
(0) calculates Kb for surfactant films, analyzing the role 
of hydrophobic tails, as a function of the tail length and 
the area per molecule, and finds a wide range of values 
for k b [Kb = with k from ref. (@)) including val- 
ues small enough, corresponding to soft interfaces, to be 
balanced by our result. Thus it is conceivable that there 
can be small tubes formed from soft membranes in water 
for which the bending forces tending to expand the ra- 
dius are compensated by the Casimir attraction, and the 
tube is stabilized by sub-leading 0(1/ R 2 ) forces. 
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